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ABSTRACT 

We develop an analytic approach to study inhomogeneous reionization on large scales 
by solving the equations of ionization balance and radiative transfer to first order in 
perturbations. Given the spatial distribution and spectrum of the ionizing sources, 
our formalism can be used to predict the large scale power spectra of fluctuations in 
the abundances of HII, HI and radiation. Our approach avoids common approxima- 
tions/assumptions in existing analytic methods - for instance, we do not assume a 
specific ionization topology from the outset; nor do we make a step-function bubble- 
like approximation to the HII distribution. Applying our formalism to sources biased 
according to the Press-Schechter prescription, we find: 1. reionization always proceeds 
"inside-out", with dense regions more highly ionized, at least on large scales; 2. on 
sufficiently large scales, HII, HI and radiation exhibits a scale independent bias relative 
to dark matter; 3. the bias is suppressed on scales comparable to or smaller than the 
mean free path of the ionizing photons; 4. if the ionizing source spectrum is sufficiently 
soft, the HII bias closely tracks the source bias for most of the reionization process but 
drops precipitously after percolation; 5. if the ionizing source spectrum is sufficiently 
hard, the HII bias drops in a more steady fashion throughout the reionization pro- 
cess. The tools developed here will be useful for interpreting future 21 cm, CMB and 
Lyman-alpha forest observations, both to learn about the reionization astrophysics 
(such as the hardness of the source spectrum and therefore the nature of the ionizing 
sources) and to possibly extract interesting cosmological information. 

Key words: cosmology: theory - intergalactic medium - diffuse radiation 



1 INTRODUCTION 

After the CMB photons decoupled at redshift about 1100, 
the intergalactic medium (IGM) remained neutral until 
the first generation of luminous sources produced ionizing 
photons. Recent measurements of the spectra of high red- 
shift quasars in the Sloan Digital Sky S urvey dBecker et al.l 
l200ll : IWhite et alj|200 j ; iFan et al.ll200fj ) and of polarization 
anisotropies in the cosmic microwave background (CMB) 
by the Wilkinson Microwave Anisotro py Probe (WMAP) 
jPage et alj l200fj| ; ISpergel et ail 120061 ) indicate that the 
IGM was reionized during redshifts z « 6—15. The 
history of cosmic reionization contains abundant informa- 
tion about the formation of the first cosmic structures, 
and can be probed in various future obse r vations, such 
as redshifted 21cm signatures (|Fieldl 1 19581 ; IScott fc Reesl 
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Il99d : iMadau et ail Il997l : IZaldarriaga et all |2004). the 
kinetic Sunyaev-Zel'dovich effect (|Sunvaev fc Zel'dovicbl 
| l980 |: lGruzinov fc Hull998l;lKnox et al.lll998l;IValageas et all 
l200ll ; ISantos et al.l 12003 : ISalvaterra et al.l I2005T ). improved 
measurements of the large -angle CMB polarizat ion fluctu 



ations (|Kaplinghat et a.1.1 120031 : IZaldarriagalll997l). the ther- 
mal state of the inte rgalactic medium l|Theuns et al.ll2002l : 



jHaiman & Spaansl 


19991: lHaiman 20021: Santos et alj|2004l; 


Rhoads & Malhotra 


20041; Haiman & Cenll2005l). Being di- 



rectly related to many observables, the distribution of the 
HI/HII regions during reionization has been studied ex- 
tensively using both semi-analytic models and hydrody- 
namic simulations (see, e.g., the recent review by Choud- 
hury & Ferrara (2006), and references therein). If the spec- 
trum of the first ionizing sources is soft (such as a nor- 
mal stellar spectrum), then the distribution of the ionized 
regions can be described by discrete HII bubbles around 
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high density regions l|Haiman fc LoebH l997). On the other 
hand, if the source spectrum is hard and extends to X- 
ray energies (such as for miniquasars) , the photons can 
more readily escape i nto the IGM, and may ionize the low- 
density regions first dMiralda-Escude et aljlioool ; I5"n]|200ll ; 
IVenkatesan et al.ll200ll b 

The reionization topology in the above two limiting 
cases are often referred to as "inside-out" and "outside-in" 
respectively. We will adopt this distinctive terminology here, 
and we will use it to describe the ionization topology on ar- 
bitrarily large scales. Existing analytic/semi-analytic mod- 
els have generally described either one or the other of these 
two limiting scenarios (i.e. by following the filling factor of 
ionized bubbles, or by assuming a uniformly rising ioniz- 
ing background). The previous methods used to implement 
such models cannot be easily generalized to derive statistics 
of the ionization topology, as a continuous function of the 
hardness of the source spectra. The goal of the present pa- 
per is to construct a model of inhomogeneous reionization in 
which we can quantify properties of the ionization topology 
for sources with different spectral hardness and clustering 
properties. The hope is to derive from first principles, rather 
than assume, the large scale ionization topology and statis- 
tics. More generally, the nature - the spectrum, spatial 
clustering, and evolution - of the ionizing source popula- 
tion is poorly known at present, and future observations (as 
mentioned above) promise to shed much light on these quan- 
tities. Parameterized models, such as the one presented here, 
will be useful for interpreting these future observations. 

In principle, hydrodynamic simulations offer an alter- 
native and more reliable way of studying the topology of 
reionization. However, it is difficult for the state-of-art sim- 
ulations to cover large scales (^ 20Mpc) while resolving 
the sources and the gas distribution on small scales. Re- 
cently, several groups have attempted to study the large 
scale properties of reionization using large scale simula- 
tions (Kohler et al.l|2005l; llliev et al.ll2005l ; IZahn et al.ll2006l ; 



iMellema et all 12009 ; llliev et al.ll2006t ). Kohler et al. (2005) 



for instance adopt a hybrid approach which combines small 
scale and large scale simulations. As will become clear below, 
our analytic treatment is in some sense similar: we use linear 
perturbation theory to address large scale questions while 
taking into account the effects of small scale clumping in 
the background evolution. Even if in the future simulations 
that simultaneously resolve the necessary small scale struc- 
ture and span over large scales were to become available, 
we would still need an analytic framework to better under- 
stand the physics of reionization. Furthermore, the ionizing 
sources have to be inserted into the simulations, and their 
properties specified, essentially by hand. In practice, it is 
likely that the intepretation of future data on the ionization 
topology will require an exploration of various parameters 
of the sources that cannot be computed ab-initio. This, in 
turn, will require a semi-analytic, computationally less ex- 
pensive model. 

In a hierarchical universe such as our own, one ex- 
pects linear perturbation theory to work on sufficiently large 
scales. The spatial fluctuations in the ionized/neutral hy- 



1 This is by no means obvious, considering the fact that reioniza- 
tion is a messy process on small scales, with large fluctuations in 



drogen and in the ionizing radiation can be related to the 
dark matter distribution via bias factors which are scale de- 
pendent in general. For a given source distribution, these 
bias factors are determined by the radiative transfer equa- 
tions and the equation of the photo-ionization balance. The 
approach followed in this paper is to solve the linearized ver- 
sions of these equations in Fourier space. Given the source 
distributions and spectra, this approach allows us to calcu- 
late the linear biases of the ionized/neutral hydrogen and 
of the ionizing radiation, and to follow their evolution. The 
price paid is that our predictions for the spatial fluctuations 
are invalid on small scales. 

Our calculation takes into account all of the relevant 
physical processes, including photo-ionization, recombina- 
tion, the diffusion of photons, the peculiar velocity of the 
baryons, and the redshifting of photons due to the expan- 
sion of the universe. It is worth noting that in the existing 
semi-analytic models, typically only photo-ionization and re- 
combination are treated in an exact mannei0. Most of the 
other processes are either missing or treated approximately. 
Our formalism accounts for all the relevant processes (al- 
beit in a perturbative manner) and allows a comprehensive 
study of the dynamical evolution of the HI/HII regions and 
the radiation field with a general source spectra. It is our 
hope that the approach taken here can be developed fur- 
ther in future work, and will ultimately make it possible to 
constrain the high redshift source properties using future 
21cm and CMB observations. Our approach is also useful 
for addressing the question of what robust cosmological in- 
formation one can obtain from future 21cm and CMB obser- 
vations (for instance, to what extent are the various linear 
biases scale independent). 

The rest of this paper is organized as follows. In i}2] we 
introduce our formalism, and present techniques for solving 
the equations. In Sj3] we discuss our assumptions about the 
ionizing source population, based on the extended Press- 
Schechter model and assuming different source spectra. In 
S|4] we present our main results, i.e. the evolution of the 
ionization fluctuations on different scales, and for different 
source spectra. In SJS] we discuss various caveats, and possi- 
ble extensions of the present paper to future work. Finally, 
in Sj6] we summarize our conclusions and the implications of 
this work. Appendix A presents an analytic solution (up to 
the numerical solution of a simple integral equation) to the 
first order radiative transfer and ionization equations. Ap- 



the form of HII bubbles. We give some plausibility arguments in 
Appendix B. Note that exactly the same issue can be raised in the 
perhaps more familiar context of large scale structure formation 
— there, it is known that linear perturbation theory works well on 
large scales, even in the late universe (such as today) when highly 
nonlinear structures exist on small scales. Why this should be so is 
not completely understood, but some plausibility arugments were 
given by Peebles (1980). We borrow his arugments and translate 
them into the language of reionization, and discuss the conditions 
(e.g. scales) under which perturbation theory is expected to work, 
in Appendix B. 

2 However, we note that Chiu et al. (2003) has a more elabo- 
rate model for the evolution of the spatially averaged ionization 
fraction, but they do not discuss the fluctuations, and/or the de- 
pendence of the fluctuations on the source spectrum. 
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pendix B contains a preliminary discussion of the validity of 
linear perturbation theory. 

Throughout this paper, we adopt a standard ACDM 
cosmological model, with the parameters fl m — 0.3, £7a = 
0.7, Q b = 0.048, h = 0.69, n = 0.95, a 8 = 0.826, 
and r = 0.088, favored by the combination of the three- 
year WMAP data and the weak lensing data of Canada- 
France-Hawaii Tel escope Legacy Survey (CFHTLS) (see 
ISpergel et aT]|2006l ). 



by Hel, rather than HI). We postpone the study of helium 
reionization to future work. Taking into account peculiar ve- 
locities, photo-ionization and recombination, the equation 
for ionization equilibrium is given by 



dn H n 
dr 



+ V • (uhiiu) 



[uh — nun , 



a B n HII 



(/// 



(3) 



rf 2 0n 7 ^\ K(fi, (f) 



a 2 W 



2 THE FORMALISM 

We begin with definitions of the physical quantities and de- 
scriptions of relevant equations in H2.ll The equations are 
solved by splitting into two pieces: a spatially averaged piece 
in H2.2l and a (first order) fluctuating piece in 



2.1 Definitions and Basic Equations 

Let us first define the relevant physical quantities: 

riHii = n H n(x,r) (1) 

riH = n H (x,r) 

n 7 = n~f(S,T,n,Cl) 
S = S(x, t, n, fi) 

where uhii and uh are the comoving number densities of 
the ionized and total hydrogen atoms (ionized+neutral) re- 
spectively. They are defined as functions of the comoving 
coordinate x and the conformal time r. Here, n 7 refers to 
the comoving photon number density per unit solid angle 
d 2 Q around the propagation direction Q. and per unit fre- 
quency parameter y,. The frequency parameter fi is defined 



H — In v — In i/q 



(2) 



where v is the photon frequency, uq = 13. 6e V/ (2nK) is the 
ionization threshold of hydrogen, and Ti is the Planck con- 
stant. The frequency parameter will turn out to be more 
convenient to use than the frequency itself. The quantity 
S/4ir is the differential ionizing emissivity, which gives the 
number of photons emitted by sources per unit comoving 
volume, per unit conformal time, per unit frequency param- 
eter fj, and per unit solid angle. 

It is useful to relate n 7 and S here to perhaps more 
familiar quantities. The proper specific intensity J of the 
ionizing radiation is related to n 1 via 

J = hn~f/a 3 

where a is the scale factor. Note that we set the speed of light 
to unity throughout this paper. The quantity S is related to 
the proper emissivity j as usually defined via 

j = TiS/{ATia A ) 

Throughout this paper, for simplicity, we will ignore 
the presence of helium atoma^. This could affect our re- 
sults somewhat for the hard-spectrum cases considered be- 
low (since most of the > 100 eV photons will be absorbed 

3 Note that we do not replace helium with hydrogen i.e. we have 
riH = Yh * Qb/mp, with Yjj = 0.76 the hydrogen mass fraction. 



where a(r) is the cosmological scale factor; u is the comoving 
velocity of the ionized hydrogen atoms; is the photo- 

ionization cross section (Osterbrock et al. 2005); and ub ~ 
2.6 x 10 _13 cm 3 s _1 is the case B recombination coefficient at 
temperature equal to 10 4 -K\ We have implicitly assumed 
electric neutrality everywhere in the universe, and we are 
ignoring the electrons that would result from the ionization 
of helium, therefore the recombination term is proportional 
to n 2 HII . The factor rcfjti, <j>) = 1 + C(exp(/x) - 1)(1 - (p a ) b 
is included to account for multiple ionizations by an X-ray 
photon through secondary ionizations by the fast photo- 
electrons. Here <f> = uhii/tih is the local ionization frac- 
tion, and we adopt the values of C = 0.3908, a — 0.4092, 
an d b = 1.7592 in the fittin g formula above (according 
to lShull fc Van SteenberSligsyh . The evolution of the ra- 
diation background is affected by the sources, the photo- 
ionization process, the diffusion of photon^ and the red- 
shifting due to the expansion of the universe, all of which 
are reflected in the follow ing radiative transfer equation (e.g. 
iGnedin fc Ostriked[l997h IT. 



Or Ofi 
S a(fi) 

(TIH - UHIljUy , 

4-7T CT(t) 



(4) 



where H(t) 



is the Hubble parameter. Our main 



task is to solve eq.[3] and eq.[4]. For our purpose, it is useful 



4 The recombination coefficient is weakly temperature depen- 
dent. At temperatures T typical of the photoionized intergalactic 
medium (where recombination is relevant) , the recombination co- 
efficient scales roughly as T -0 7 . The temperature is in turn re- 
lated to gas density to some power, with the power index ranging 
from about 0.0 to 0.6 l|Hui fc G ncdin 1997). Overall, the spatial 
fluctuation of the recombination coefficient in ionized regions due 
to the fluctuation in density (and therefore temperature) is rather 
weak, and is ignored here. 

5 We caution that the fitting formula only works well for high 
energy photons ( > lkeV). 

6 Note that the word "diffusion" refers to the term f2 ■ Vn 7 in 
eq. [4], and does not imply scattering. This is a very loose usage 
because this is not the usual term in the diffusion equation. 

7 This equation ignores recombination radiation, gravitational 
redshift and Doppler shift by peculiar motion. Gravitational red- 
shift is a negligible effect except on scales comparable to the hori- 
zon. Doppler shift by peculiar motion does not contribute to first 
order in perturbations (after averaging over the photon directions, 
which is what we will do eventually); its contribution to second 
order is negligible compared to other existing second order terms. 
The importance of recombination radiation is diminished to some 
extent by the cosmoloigcal redshift. 
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to rewrite the definitions in eq. [T] in terms of the spatial 
averages and the perturbations: 

n H u = n H fHii(T)[l + Shii(S, t)) (5) 
= n H [fan (t) + A H ii (x, r)] 

nn = uh[1 + S(x,t)] 

n 7 = UH f-y(T, fj,)[l + 5~/(x,T, fJ.,0,)) 

= n H [fj(T,fi) + A 1 (x,T,^i,Q,)] 

S = n H f s (T,fi)[l + 6 s (x,t,h, n)] 
= n H [fs{r,ii) + A s (x,t,/i,Q.)] 

where fmi and / 7 are the mean number densities of ion- 
ized hydrogen atoms and photons respectively, and / s is the 
mean source emissivity - all normalized by the mean co- 
moving total (i.e. ionized plus neutral) hydrogen number 
density uh, which is a constant in time. S, Shii, 5 7 and S s 
are the corresponding overdensities. Ahii, A 7 and A s are 
introduced because they greatly simplify the following calcu- 
lations. We will assume below that the (total) hydrogen fluc- 
tuations faithfully trace the dark matter fluctuations (with 
no bias), which is .justified on scales well above the Jeans 
scale (e.g lGnedin fc Huilll998l ). Note that according to the 
above definitions the neutral hydrogen density is given by 

riHi = n H [(l - fmi) + 8 - Ami] 



and 



riHi — riHi 
uhi 



8-Ai 



fmi 



(6) 



2.2 The Spatial Averages 

2.2.1 The Exact Solutions 

Taking the spatial (and angular) averages of eq. [3] and eq. [4] , 
we find: 



Or 



= Ml-f H ii) I dn^(K)f y Ci x H 



■,(D 



dr 



OLEUM 2 n 
o JHII^HII 



— 5-(l - fmiJfaC H 



(8) 



where 0$ and Cfi are the clumping factors for photo- 
ionization, defined as: 



a 



(1) _ (UHin-fK) 



a 



"< H {nm){nJ{ K ) 

(2) _ (n gJ n 7 ) 
■yH 



0) 



(nm)(n-f) 

Chii = {n H ii) I '{nun) 2 is the clumping factor for recom- 
bination. To solve for fan and / 7 , one needs to know the 
evolutions of C^ H , C^ H , and Chii- Since the clumping fac- 
tors are likely dominated by non-linear density variations 
on small scales (e.g. Haiman, Abel & Madau 2001) their 



values cannot be computed reliably in our present frame- 
work. We assume, for simplicity, that the clumping factors 
are known from small scale hydrodynamic simulations. Note 
that recent large scale radiative transf er simulations effec - 
tively adopt the same assumption (e.g. iKohler et al] [2005). 

Instead of r, let us use lo — In a(r) as the time vari- 
able. To further simplify the notation, we define <j(t, /i) = 
a(n)fiii/{Ha 3 ) and Qs(r) = olbtih /(Ha 3 ), which can be 
interpreted as follows: a is the probability that a photon 
of energy p, experiences a direct photo-ionization in Hubble 
time if the universe is neutral; &b is the average number of 
recombinations a proton would experience in a Hubble time 
if the universe is completely ionized. Eq.[7] and eq.[8] can 
then be rewritten as: 



dfmi 



4tt(1 - f H ii) 

OBfHIlCHII 



d^<«)/7C. (1) 



-II 



dfy _ fs £>jj_ 
duj 4nHa dfi 



(10) 



(11) 



- (1 - fmiYof-,C^ H 

Eq.[lO] and eq.pTJ can be solved iteratively. To do so, we 
need to use a trick which also appears in ^2.31 First, let us 
change the variables from (u>,fj,) to (u,v) which are defined 
as: 

1 



u = -(u + p.) 



v = 



(12) 



Therefore 
d d d 
dv dcj dp 

Eq.pT] can be transformed into: 

9fj fs /, f \ ~ c n (2) 

or in a simpler form as: 
Of 

— = q(u, v) - p(u, «)/ 7 



(13) 



(14) 



(15) 



where q = / s /(47riifa) and p = (1 — fni^oC^. Eq.[l5] is a 
standard first order differential equation. Its solution reads: 

dv'q(u, v') exp[— / dv"p(u, v")] (16) 

■ oo J v' 

or in terms of the variables (oJ,p) as: 

dio'q(uj', p + lo — lo') (17) 

- oo 

x exp[— / dui"p(uj" , p + lo — lo")] 

The above solution assumes that there are no ionizing pho- 
tons at arbitrarily early times and/or arbitrarily high ener- 
gies. 

Using eq.pTj, we can obtain / 7 (o;,/i) for an initial 
fHii{to)- f y can then be used to calculate a new fun using 
eq.JTO). In practice, we find that if this procedure is repeated, 
/ 7 and fun converge to their true values after around ten 
iterations. 
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2.2.2 Useful Approximations 

Eq.[l7] presents a useful picture of physics: it shows that at 
the early stage of reionization (when the IGM is not highly 
ionized), the low energy photons (but above 13.6eV) have a 
short memory of their past history because they are quickly 
absorbed by neutral hydrogen. This can be seen from the 
largeness of the exponent in eq.JTT] when the cross section 
is large. In contrast, the hard photons retain a long mem- 
ory. Another way of saying the same thing is that the soft 
photons have a short mean free path compared to the hard 
photons. This fact leads to a much simpler way of solving 
eq.[TD] and eq.QT] when the average ionized fraction fnn 
is not close to one. First of all, we notice that in the low 
energy limit, eq.JTT] reduces to a simple form: 



u 



(18) 



4tt(1 - f H ii)HaZC™ 



Secondly, we assume there is a critical frequency parameter 
fl c , above which the photons hardly ionize any neutral hy- 
drogen during reionization, and below which eq. [18] is valid. 
Therefore, eg. [10] can be rewritten as: 



dcu 



aBCmifHii 



(19) 



Eq. [18] essentially describes an emission-absorption equilib- 
rium, i.e. all emitted photons with < fi < p c are consumed 
by ionization. It should be noted that if the source spectrum 
is hard (i.e. high percentage of photo-ionization is caused by 
high energy photons.), Eq. [18] and [19] are not a good ap- 
proximation anymore. This will be discussed further in £|2.3I 
and |4] Note that throughout this paper, we compute the 
exact solutions rather than employ the approximations out- 
lined above, though we will compare the two in U 



2.3 The Linear Perturbations 

2.3.1 The Exact Solutions 

Let us denote the Fourier transforms of S, A s , Ahii and A 7 
as S, A s , Ahii and A 7 respectively. The Fourier transforms 
of eq.[3] and eq.[4j to the first order are: 



dAi 



dtu 

dA 7 _ 
dcu 

where 
F = 



= G5 
<9A 7 



FAhii + 



dfj,{K) 



d SIA-/B 



MA 7 + NA a + R{Ahii - 5) 



+ 



G = 



26>BfHII 
!• oo 

4-7T / dfiaf. 
Jo 

dlnD „ 



(20) 
(21) 

(22) 



(«> - (1 - fHIl) g 



4>~fHII 



du> 



JHII 



8 The question of whether /when retaining only first order per- 
turbations is justified is discussed below and in Appendix B. 



+ 


f 00 

4tt / d/xaf- 




J 


S3 — 


(1 - Jhii ja- 


M = 


il - Smi)a - 


TV = 


(AixHa)- 1 


R = 





(K) - (1 - f H n)fB 



ik ■ Q 



c)k 



4>=!hii 



Ha 



and k is the wave vector. For simplicity, we do not show 
the k dependence explicitly for the Fourier modes. In deriv- 
ing eq.[20], we have used the fact that to the first order, on 
large scales, the dark matter overdensity grows linearly with 
a growth factor defined by D(t), and the peculiar velocity 
is proportional to the gradient of the gravitational poten- 
tial. The great advantage of the Fourier transformation is 
that it allows us to solve eq. [21] using tricks similar to those 
introduced in the previous section, i.e. : 



A 7 (w,^,!T2) = / du'{N(u')A s {u' + -w',fi) (23) 

J — OO 

+R(lu', li + lu- u')[A H ii(lu') - S(lu')]} 



x exp[- 



dio" M(lu" , n + lu — uj" , fl)] 



The first line above describes the contributions to fluctua- 
tions in radiation from fluctuations in the spatial distribu- 
tion of the ionizing sources; the second line describes the 
contributions from fluctuations in absorbers; the third line 
accounts for the modulation by optical depth, i.e. it tells us 
the distance to which one needs to integrate. By integrating 
eq.[23] over Cl, one can obtain the monopole perturbation of 
the radiation field, which is what is needed in eq.|20]: 



d 2 f2A 7 (tj, /x, D.) 
= 4tt / du'{N((j')A a (u>',ij, + L>-u') 

J — OO 

+ fl(w',p + W-w')[&HX/(w') -*("')]} 

x exp[— J duj"B(uj", ji + uj — lu")] 
sm\P(u,u>')k] 



(24) 



P(uj,io')k 



where 



P(u,u') = 



dui" 



H(ui")a(uj") 



(25) 



In writing down the above expression, we have assumed 
that the dominant contribution is from the monopole of A s 
(note that we have dropped the argument Q, to signify the 
fact that this is the monopole component). Note that we 
are not assuming A 7 has no higher multipoles; rather, we 
are assuming that, as far as the source contribution is con- 
cerned, the monopole of A 7 is dominated by the monopole 
of the source emissivity. This assumption can be motivated 
in two different ways. First, in the low k limit (which is the 
regime where perturbation theory probably works best), one 
can show that the monopole dominates. Second, on large 

9 Throughout this paper, we approximate (k((j,, cf>j) by «(/i, {4>))- 
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scales, since one is averaging over many sources, the result- 
ing smoothed emissivity is probably close to isotropic even 
if the individual sources are not. 

Eq. [24], together with eq. [20] , allows for an iterative so- 
lution for Ahii and A 7 given the source distribution A s and 
the dark matter overdensity 8. In Appendix A, we present a 
further improved scheme of finding the numerical solution. 

2.3.2 Useful Approximations 

Similar to what we have done in £|2-2I when the universe 
is not close to being fully ionized, eq . [24] for soft photons 
can be greatly simplified. First, we notice that the factor 
sin[P(uj, u)')k]/[P(ijj, ui')k] in eq.[24] should be very close to 
unity for soft photons. There are two reasons for this: one 
is because the wave number k of interest is small; the other 
reason is that the large exponent B limits the magnitude of 
co — u)' to be very small in the integration. Taking these into 
account, eq.[24] becomes much simpler: 



(26) 



4tt 



{N(u>)A s (u> >f j.) + Biu,p)\h B xi{pi) - S(u)]} 



The above expression is the analog of eq. [18] . representing 
essentially emission-absorption equilibrium. Using the same 
critical frequency parameter fi c to isolate the contributions 
from the soft photons as in eq.[l9], and assuming eq.[26] is 
correct for such soft photons, we obtain: 



dA f 



/dln_D _ \ 



du) 

{2&BfHII 



- V) fmiS 
Y)A HII 



(27) 



+ 



— J dfi(K)A s (uj,n) 



where 



Y = 4tt(1 - f h 



(28) 



4>=1hii 



Again, eq.[27] is a good approximation when the source spec- 
trum is not too hard. In ^4] we quantify this statement and 
discuss the choice of fj, c with realistic examples. We reiter- 
ate that all of our conclusions in this paper are based on 
the exact solutions rather than the approximations outlined 
above, though we do compare the two in §4] 

At this point, the reader might wonder: since reioniza- 
tion is likely a complicated process with large fluctuations on 
small scales, could these fluctuations invalidate linear per- 
turbation theory on large scales? In other words, how should 
one justify the use of perturbation theory on large scales? 
This is actually a difficult question. The same question arises 
in the context of large scale structure formation: how do we 
justify the use of linear perturbation theory on large scales 
today, knowing that there are highly nonlinear structures 
on small scales, such as gala xies, clusters and so on? Some 
plausibility arguments exist (|Peebleslll980h . and we will dis- 
cuss the analogs of these arguments for reionization in Ap- 
pendix B. However, the only rigorous method of valdiation 
that we know of is to appeal to numerical simulations. We 
will perform comparisons of our calculations with numerical 
radiative transfer simulations in another paper. 



At least one difference between reionization and large 
scale structure is, however, worth emphasizing. In the case 
of reionization, the spatial averages are undoubtedly af- 
fected by small scale clumping - hence the need to intro- 
duce clumping factors in eq. [7] and |5j. In some sense, our 
treatment here i s quite similar to so me of the recent numeri- 
cal simulations l|Kohler et al.ll2005l ) which incorporate small 
scale clumping by hand while focusing on the large scale 
fluctuations. In the case of large scale structure, the spa- 
tially averaged equations are those that govern the global 
expansion of the universe i.e. the Friedmann equation and 
energy-momentum conservation. In that case, small scale 
clumping appears not to a ffect significantly t he evolution 
of the spatial av erages (e.g. iHui fc Selial3 ri996: see however 
iKolb et al.ll2005l for a different view). 



3 MODELING THE IONIZING SOURCE 
POPULATION 

Solving the equations described in the previous section re- 
quires specifying the source properties, including the emis- 
sivity as a function of redshift, and the spatial distribu- 
tion and spectrum of the ionizing source s. In this paper, we 
use t h e extended Press- S chechter theory llPress fc Schechterl 
1 19741 : iBond et al l Il99ll : lLacev fc Coll Il993l ) to model the 
abundance and spatial distribution of the ionizing sources. 
It is important to emphasize that this is for illustration only 
- our formalism as laid out in the previous section is cer- 
tainly not wedded to this particular model of the ionizing 
source population. 



3.1 Dark Matter Halo Abundance 

The minimum halo that can host luminous sources should 
have a virial temperature abov e 10 4 K to allow effici ent hy- 
drogen line cooling (see, e.g. iMesinger et al.l l2006). This 
leads to a minimum halo mass given by: 



0.3 



II 

- 1 / 2 / h \ 1 I „... , A - 2 



M mm « l.3xlO 7 M (^) (i±£)" (29) 



\0.7J V 1.22 J 



where jj, mo i is the mean molecular weight, which is chosen 
to be 0.6 (appropriate for the ionized gas in halos with a 
virial temperature above 10 4 K) in the following calculations. 
According to the extended Press-Schechter Model, on mass 
scale m, the fraction of mass collapsed in halos with masses 
larger than Mmin is: 



_C"(f,r)=erfc 



8 C — 5 m (x, t) 



(30) 



where S c is the critical overdensity in the spherical collapse 
model; dm and axe the overdensity and the variance of 
the density fluctuations on mass scale m respectively; is 
the density variance corresponding to the mass scale Mmin. 
On large scales, a m is much smaller than o"^ in and therefore 
neglected in the following calculations. 

Assuming on average each hydrogen atom in the col- 
lapsed objects emits y(fj,) ionizing photons per unit fre- 
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quency parameter /j, the emissivity function smoothed over 
scale m is given by: 

d 

S m (x,T,fi) =l{n)n H -Q^ [fm ll (x,r)(l + 8m(x,T))] (31) 

The factor 1 + <5 m takes into account the mass overden- 
sity. This equation also implicitly assumes that the collapsed 
objects produce their photon output on a time-scale much 
shorter than f co "/[df co "/dt] (which is justified for short- 
lived, massive stars, or efficiently accreting black holes). If 
we choose the scale m to be very large (e.g. the horizon 
size), 5 m can be neglected, and from eq. [31] we obtain the 
spatially averaged emissivity function defined in Sj3]as: 



(32) 



f 51 ^ 










V V 2(T Ln J 



By Taylor expanding eq.[3T] around 8 m — and then per- 
forming a Fourier transformation, we obtain the spatial fluc- 
tuation of the emissivity as ™ 



A 3 (fc,r )M ) = -y(^^-[R(r)S(k,r)] 



dr 
8(k,r) d 
D(r) 8t 



(33) 



[R(t)D(t)] 



where 



R(t) = erfc 



5 C 



2 

min 



+ 



■ exp 



1°l 



(34) 



and D(t) is the linear growth factor. 

In our model, we do not take into a ccount the 
feedback effects discussed in recen t papers dCole et all 



200d: lot] l200ll: iBenson et all Eioi lOh fc Haimanl 
Diikstra et al.ll2004l; [Benson et al.ll2005l ; iKramer et all 



2003 



2006 



Mesinger et al. 20061 . also see Haiman & Holder 2003 for 



a general reference for feedback effects during reioniza- 
tion). For example, the ionized regions usually have rela- 
tively higher temperat ures due to photo-he ating, leading to 
a larger filtering scale l|Gnedin fc Hui|[l998h that suppresses 
the formation of small scale structures. This reduces the 
clumping factor and hence the recombination rate, and also 
results in an anti-correlation between the source overdensity 
and the local ionized fraction. In principle, such a feedback 
effect can be included in the linear perturbation calculation 
by inserting into eq. [33] a term that is proportional to Ann- 
We leave a more careful study of such effects to the future. 

In our calculation, the normalization of the source emis- 
sivity is chosen to give a mean optical depth of Thomson 
scattering equal to 0.088, which is suggested by the three- 
year WMAP data. The remaining freedom is to choose the 
sourc e spectrum , which depends on the type of sources (see, 
e.g. . IOhll200ll ). Rather than exhausting spectra of general 
shapes, we focus on spectra of power-law forms with different 
spectral indices, which cover the possible range of effective 
spectra of stars and miniquasars. A power-law form in the 
frequency v is an exponential form in the frequency param- 
eter fi, i.e. : 



7(A*)d|i = exp[(/3 + l)fj]dn 



(35) 



where j3 is the spectral index, Cp is a normalization factor, 
and C, is the total number of ionizing photons generated per 
baryon in stars, and managing to escape into the IGM. For 
example, if 10 percent of the baryons turn into stars with 
a normal Salpeter mass function, and 10 percent of their 
ionizing radiation escapes, then £ = 40. The spectrum is 
smoothly cut off at n = 10 (~ 300 keV) to allow for a 
proper normalization even when /3 ^ — 1. The cutoff does 
not introduce any artificial effects because the mean free 
path of photons at this energy greatly exceeds the Hubble 
distance. 

It is worth noting that in a more complicated scenario, 
the shape of the source spectrum may vary with redshift. 
For example, if quasars are the dominant ionizing sources 
at high redshifts (z ~ 15), and followed by stars at z ~ 6, 
the hardness of the source spectrum varies with time. Such 
a case will be studied in a future paper. 



4 RESULTS 

As we have discussed in H2.2I the spatial averages of the 
ionized fraction and the radiation intensity can be calcu- 
lated by solving eq.|10| and eq.[8] once the clumping factors 
C^ H , C^j, and Chii are provided. According to the ex- 
isting hydrodynamic simulations, Chii is around order of 
ten during reionization for a UV dominated source spec- 
trum, but can be less than one for a harder source spec- 
trurrrH. The values of CjQ and CfQ vary with both the 
redshift and the ph oton energy, but are rarely far from unity 
|Kohler et al.ll2005T ) . For simplicity, in this paper, we choose 
Chu = 10 for (3 = -3 or (3 = -2, Chii = 1 for j3 = -1, and 



C^ H = G K * H = 1 for all cases. We caution that the precise 
values of the clumping factors are rather uncertain. The im- 
portant point to keep in mind is that the clumping factors 
show up explicitly only in the equations for the spatial aver- 
ages but not for the fluctuations. In other words, they affect 
the linear perturbations only indirectly through their effects 
on the background (i.e. fnn and / 7 ). Once the background 
is specified, our predictions for the linear perturbations are 
quite robust. This is discussed further in SJS] below. 

To illustrate the results, we choose three types of source 
spectra: (( = 82, fi = -3), (C = 75, /3 = -2) and (C = 
50, f3 = —1). We show the evolution of the spatially aver- 
aged ionization fraction fnn in Fig. [T] the HII bias (i.e. 
5 H ii(k)/5(k)) at k = 0.01 Mpc -1 in Fig.H and the HII bias 
as a function of the scale k at redshift z = 20, 13, 10, 9 in Fig. 
[3] From Fig. [2] we see that in the case of a soft source spec- 
trum (i.e. f3 = —3 and f3 = —2), the bias of the HII regions 
remains at a high value during most time of reionization, 
and quickly drops to one when the HII regions percolate the 
IGM. This can be understood in the usual bubble picture, 
in which the HII regions are confined within the HII bubbles 
due to the short mean free path of the photons. Therefore, 



10 Note that the wave mode of the Fourier transformation al- 
ready indicates a smoothing scale, therefore the scale index m is 
dropped. 



11 Note that the last statement is probably not true at the end 
of reionization when even high density regions are ionized and 
Chii 3* 1. For simplicity, we ignore the time dependence of the 
clumping factors in this paper. 
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Figure 1. The solid curves show the evolution of the ionized 
fraction fun\ the dotted curve is the approximate solution for 
fun from eg. 1191 . The three panels assume different power law 
slopes (/9) for the source spectrum, as labeled. 



Figure 2. The solid curves show the bias factors of the HII re- 
gions at k = 0.01 Mpc -1 ; the dotted curves show the approximate 
solution from eg. 1271 . &hii an d S represent the Fourier transforms 
of the overdensities Suil an d S, defined in eq. [5], 



the HII regions do not merge on large scales until the mean 
ionization fraction of the IGM reaches a very high level. In 
contrast, for a harder source spectrum, the bias of the HII 
regions decreases in a more continuous fashion. Interestingly, 
this is not only because hard photons have a much longer 
mean free path, but also due to the fact that secondary ion- 
izations are more intense in less ionized regions. 

The HII bias is intimately related to the issue of inside- 
out versus outside-in ionization. A reasonable definition of 
inside-out ionization is that (SSx) > 0, where 8x is the 
fluctuation in the ionized fraction. A positive correlation 
means higher density regions exhibit a higher ionized frac- 
tion. Conversely, a negative correlation can be taken as 
the definition of outside-in ionization. To the lowest order, 
(SSx) = {S(Shii — S)). Therefore, the sign of (SSx) is deter- 
mined by whether the HII bias is greater or less than unity. 
It is worth noting that whether reionization is inside-out or 
outside-in can be a scale dependent question: the sign of 
(SSx) could depend on the scale of interest (think of 8 and 
5x as quantities smoothed on some scale, or think of their 
Fourier counterparts). 

An important feature we learn from Fig. [5] is that the 
large scale bias of the HII regions is larger than unity in all 
three cases, which therefore means the high density regions 
tend to be more ionized than the low density regions. This 
fact shows that at least on large scales, the topology of the 
HII regions is inside-out, even for a hard source spectrum. 
This seems to hold on all scales that we have examined (see 
Fig. [3]l , although one must note that our perturbative ap- 
proach is expected to break down on sufficiently small scales. 

By choosing fj, c = 3.75 (E v ~ 580eV), we find that 
eq.[l9] and eq.[27] are good approximations for fun and 
5hii on large scales {e.g. k <, 0.01 Mpc -1 ) during most of 



t<0 



["O 



m — 



lO — 
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0=-l 
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Figure 3. The HII bias as a function of scale. The long dashed, 
short dashed, dotted, and solid curves are for redshift z = 
20, 13, 10, and 9 respectively. The arrows indicate the mean free 
path of the typical ionizing photon in each case at z = 9 (see 
eg. 1371 ). We caution that our linear theory becomes inaccurate 
on small scales (fc > 0.1 Mpc -1 ). 
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the reionization process for the /3 — —3 and = —2 cases, as 
shown in Fig.[T]and Fig. [2] The reason for choosing this value 
of jj,c is not only that it provides a good fit for both fnu 
and Shu, but also because 580eV corresponds to a critical 
energy threshold, below which photons are significantly ab- 
sorbed by neutral hydrogen, as shown in Fig. ijPl . In other 
words, eq. [18] works well for photons with energies below 
580eV during the early stages of reionization. For smaller 
scales, the agreement becomes worse because the photon 
propagation suppresses the small scale fluctuations, thus re- 
ducing the amplitude of the (HII to dark matter/baryon) 
bias. This is shown in Fig. [3] 

A useful way to see why the bias of HII decays faster for 
a harder source spectrum is to calculate the mean free path 
of the ionizing photons. For a specific frequency parameter 
fj., the comoving mean free path is simply a 2 /[uh cr(/i)(l — 
fan)]- The average mean free path can be defined as: 

_ \ n H (l - fan) Iq° MM ] ~* , . 

" [ * xrrrJ (36) 

Using cq.[TS] for three different source spectral indices, we 
find: 




These numbers basically answer the question of why in Fig. 
[2] the HII bias in the /3 = — 3 or /3 = - 2 case does not 
drop until the ionization fraction is very close to one. This 
is because in both cases the wave length of k — 0.01 Mpc - 
is larger than the average photon mean free path during 
most time of reionization. Radiative transfer is unable to 
suppress perturbations on scales larger than the mean free 
path, and the HII bias more or less tracks the source bias. 
In contrast, HII perturbations on smaller scales are more 
easily suppressed by the radiation, which can be seen in 
Fig. [3] as the suppression in Shii/S for k ^ 2tt/X. However, 
we caution that our linear theory becomes inaccurate on 
small scales. We also find that other ways of estimating the 
average photon mean free path only change the results by a 
factor of a few, meaning that eq.|36| provides a robust order 
of magnitude estimate. For example, if secondary ionization 
is included in eq. [36], the average mean free path in all thre 
cases is increased by no more than a factor of two. 

In Fig. 2] we show the redshift dependence of the fluctu- 
ation of the neutral fraction with respect to the dark matter 
(i.e. (1 — fHii)Sm /S) at k = 0.01 Mpc -1 , which is propor- 
tional to the 21cm signal f rom the neutral hydrogen (see 
e.g. IZaldarriaga et al.ll2004l : IZahn et al.ll2006l '). At the early 
epoch of reionization, this quantity is close to unity since 
most of the IGM is still neutral. The non-monotonic behav- 
ior at lower redshifts may suggest a best window for detect- 
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Figure 4. The redshift dependence of (1 — fHll)&Hl/& at ft = 
0.01 Mpc" 1 . 

ing the 21cm emission. Such a feature is more pronounced 
when the source spectrum is soft. 

The bias of the radiation monopole (J d 2 fW 7 /(47r<)), 
which is simply called 5^/5 in the figures) on scales of 
k = 0.01 Mpc -1 and k = 0.1 Mpc -1 is shown in Fig. [6] as 
a function of the photon energy for different redshifts, and 
in Fig. [7] as a function of redshift for different photon ener- 
gies 13 !. Again, one can see a clear difference between the soft 
and the hard photons, which are divided by the critical line 
at E v ~ 580eV. The bias of the soft photons remains high 
until the HII percolation, following closely the bias of HII 
itself, at least during the early stages of reionization. On the 
other hand, the high energy photons diffuse relatively freely 
into the space, leading to an ever decreasing bias. 

An interesting feature in Fig. [7] deserves a brief dis- 
cussion: the bias of the radiation field (especially for soft 
photons) shoots up quickly right before percolation. In Fig. 
[2] we notice that for a soft source spectrum, the bias of 
HII also rises before the time of percolation (but less dra- 
matically). This appears counter-intuitive: naively, one ex- 
pects that both the HII abundance and the radiation roughly 
trace the ionizing sources, whose bias in our version of the 
extended Press- Schechter model always decays with time 
(shown in Fig. [8}. While this intuition is correct in the early 
stages of reionization, the situation starts to change when 
the mean ionization fraction becomes significant. This is be- 
cause a high density region generates more photons than it 
can consume (by ionization when the ionized fraction is al- 
ready high), and this leads to run-away, causing the abrupt 
rise of the bias of the radiation field (with a corresponding, 



12 Note that the y-axis in Fig. [5] is in logarithmic scale. We do 
not show labels on the y-axis because they are not important for 
our purpose. 



13 Unlike the HII bias, the radiation bias is probably unobserv- 
able. We plot it here simply for a better understanding of the 
physics of reionization. 
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but less dramatic, rise of the HII bias). After percolation, 
the radiation bias of course drops precipitously because the 
photons are free to diffuse to large distances. This effect is 
less pronounced for hard photons because percolation for 
hard photons is a more gradual process to begin with. 

Another interesting feature in Fig. [7] is that the bias of 
the soft photons exhibits damped oscillations after reioniza- 
tion is complete. This behavior can be traced back to the 
oscillatory kernel eq.[21]. In the case of j3 = —3, the bias of 
the soft photons right before the percolation of HII bubbles 
rises up to a high value, meaning that high density regions 
contain many more photons than low density regions. Such 
a difference between the high and the low density regions 
quickly decays away when the HII bubbles merge, i. e. when 
the soft photons are "released" and can freely travel to the 
neighbouring low density regions. The relaxation of this pro- 
cess leads directly to the oscillations we see in Fig. [7] The 
oscillations are most pronounced when the wavelength is 
small (high k) and when the post-percolation drop in bias 
is most abrupt (i.e. when the photons are soft). Note that 
these oscillations can be washed out by the stochastic bias of 
the ionizing sources, which is not considered in this paper. 

Finally, we find that the biases of both HII/HI and 
the radiation field remain scale invariant on large scales 
(k ^ 0.01 Mpc -1 ). This is not only because the source bias is 
scale invariant in our model, but also due to the fact that the 
diffusion of photons is negligible when the scale of interest 
is much larger than the mean free path. Using eq. [M] again, 
we can see that at small k the oscillatory kernel remains 
constant, therefore the k dependence is essentially removed. 
The scale invariant nature of the large scale HII/HI bias can 
be very useful, because it makes it possible to measure the 
shape of the primordial mass power spectrum using the HI 
power spectrum. The evolution of the HI fluctuation shown 
in Fig. [4] may be useful for this purpose in future 21cm ob- 
servations. 



5 DISCUSSION 

In this section, we briefly discuss several interesting issues 
related to our results above. First, recall that the clumping 
factors Cijj and Chii are set by hand in our calculations. 
In particular, we have neglected their dependence on red- 
shift and frequency (the latter is relevant for C&g )• The 
precise dependence is uncertain, and is only p artially ad- 
dressed by the most rece nt simulations (see e.g. Illiev et all 
120051 ; iKohler et al]|2005f ). However, as emphasized before, 
the clumping factors appear explicitly only in the evolution 
equations for the spatially averaged quantities (eq. [7] & [8]), 
and not in the equations for the first order perturbations (eq. 
|20| & [21] ) . The clumping factors only influence the first or- 
der perturbations indirectly through their influence on the 
background quantities fun and / 7 . While quantitative de- 
tails regarding the perturbations do depend on the precise 
values of the clumping factors, the overall qualitative behav- 
ior, such as the evolution of the various biases, remain quite 
robust. 

Another issue concerns the topology of the HII distribu- 
tion. It is interesting to ask if the topology is still inside-out if 
the source spectrum is made up of only hard photons. This, 
for example, is relevant to the scenario proposed by Ricotti 
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Figure 5. The spectrum of the ionizing background. The vertical 
line refers to E 7 = 580 eV. The solid, dotted and dashed curves 
are for z = 18.7, z = 13.6 and z = 9.8 respectively. The bended 
curves are the exact solutions from eq . 1111 . the straight lines are 
the approximated solutions from eq.|18|. 
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Figure 6. The bias in the ionizing background photon density as 
a function of photon energy. The vertical long dashed line refers 
to E~i = 580 eV. The solid, dotted and dashed curves are for 
z = 18.7, z = 13.6 and z = 9.8 respectively. The left column is 
for k = 0.01 Mpc -1 , and the right column is for k = 0.1 Mpc -1 . 
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Figure 7. The bias in the ionizing background photon density as 
a function of redshift. The solid, dotted and dashed curves are for 
E 1 = 170 eV, B 7 = 580 eV and _B 7 = 25keV respectively. The 
left column is for k = 0.01 Mpc -1 , and the right column is for 
k = O.lMpc -1 . 

& Ostriker (2004), in which reionization at high redshift is 
caused by highly obscured miniquasars. To answer this ques- 
tion, we consider three cases, in which we remove the soft 
photons by mutilating the = — 1 spectrum and erasing 
all photons below 170eV, 270eV", 450eV, respectively (the 
emissivity factor £ is raised to 80, 100, and 150, respectively 
to yield the same optical depth). In Fig. O we see that the 
bias of HII always exceeds unity, meaning the topology is 
not outside-in, but the HII bias does decrease with increas- 
ing hardness of the ionizing photons. One can also try to 
reduce the HII bias by lowering the source bias. In Fig. 1101 
we repeat the calculation for the case of (£ = 82, = —3) 
except that the source bias is artificially set to unity, which is 
probably the lowest value one can expect, in any scenario in 
which the ionizing sources populate collapsed halos. We find 
that the bias of HII is again always larger than oncFl. For 
= — 2 or /3 = — 1, the evolution of the HII bias becomes 
more featureless and very close to unity at all redshift f^l 
Therefore, our conclusion about the topology of the HII dis- 
tribution is robust 

14 Note that the bump at z ~ 10 is due to the run-away effect 
introduced in [|4] 

15 In the case of j3 = — 1 and a source bias of unity, the HII bias 
right before the end of reionization does go slightly below unity. 
This means that the ionization topology is marginally outside-in. 

16 We caution that our conclusion is based on the linear pertur- 
bation calculation, which does not take into account the spatial 
fluctuations of the second or higher order terms in the equations. 
The non-linear opacity and recombination fluctuations are im- 
portant in determining the shapes and the topology of the HI I 
regions on small scales (< 1 Mpc) (see, e.g. , ICiardi et al.|[200ll) . 
Whether or not the non-linear terms can affect the large scale 
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Figure 8. The evolution of the source bias <5 s /<5(= A s /(/ s <5)) 
using eg. 1321 and eg. 1331 . 



At the end of reionization, our conclusion about the 
ionization topology may seem counter-intuitive, especially 
for X-ray reionizations. Assuming a homogeneous radiation 
background (which is the limiting case for very hard spec- 
tra of the X-ray reionization scenarios considered here), one 
finds that in the limit of a low neutral fraction, s;hi = 
nm/nH, ionization equilibrium implies xbi oc tih, i.e. more 
overdense regions are less ionized. Therefore, the ionization 
topology should be "outside-in" rather than inside-out, in 
the limit of late times and a very hard spectrum. Our cal- 
culations above do not reveal this limiting case, however, 
for the following reasons: (i) our sources do not have ar- 
bitrarily hard spectra, (ii) we include a clustering of the 
ionizing sources, which, together with the rapid evolution of 
the emissivity, causes the background radiation to be per- 
sistently non-uniform. We have, however, checked that our 
code can reproduce the outside-in limiting case when we set 
the source bias to be zero (i.e. at late times, our code yields 
the expected behavior 5bii/8 ~ (1 — xbi) < 1, which fol- 
lows from xhi oc tih in the limit of a small mean neutral 
fraction xhi and small fluctuations). 

We have made at least two important simplifications 
in our calculations. One is ignoring helium. The other is to 
assume the source bias is deterministic^ It is in principle 



ionization topology is still an open question (although the argu- 
ments in Appendix B suggest that the high order effects can be 
neglected on sufficiently large scales). 

17 In other words, we have assumed the source overdensity is 
linearly proportional to the local mass overdensity. The relation 
between the two is expected to be more complex in general: it 
could be stochastic and the sources certainly have Poisson fluctu- 
ations. In essence, we assume in this paper the stochasticity and 
Poisson fluctuations are negligible on the scales of interest. 
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Figure 9. The bias of HII as a function of mean ionized fraction. 
The solid, dotted, and dashed curves correspond to energy cut- 
offs at 170eV , 270eV, and 450eV respectively, i.e. photons below 
this energy were erased for the j3 = — 1 source spectrum. 

straightforward to relax these assumptions/restrictions. We 
hope to do so in a future paper. 

In practice, we still know little about the source prop- 
erties at redshift above six. This work can serve as a link 
between the high redshift sources and the large scale dis- 
tributions of the HII regions and the radiation fields. For 
example, as we have shown, the bias of the source distribu- 
tions is similar to that of the HII regions during the early 
period of reionization; the evolution history of the bias of 
the HII regions can be used to constrain the hardness of the 
source spectrum. It allows us to constrain the source prop- 
erties with the upcoming 21cm observations and the kinetic 
SZ effect from the CMB. Furthermore, as we have found in 
our calculation that the distribution of HII traces the dark 
matter's on large scales with a scale independent bias, these 
observations may also be used to measure the shape of the 
linear matter power spectrum. A more careful discussion of 
these issues will appear in another paper. 



6 CONCLUSIONS 

We have developed a perturbation theory of cosmic reion- 
ization by solving the linearized radiative transfer equation 
and the equation of ionization balance in Fourier space. The 
formalism can be used to predict the large scale fluctuations 
of the HII regions and the radiation fields for a given spatial 
distribution and spectrum of the ionizing sources. The nu- 
merical solutions are straightforward to obtain. In the case 
of UV dominated source spectra, we have found an approx- 
imate analytic solution which works accurately in the early 
stages of reionization. 

To illustrate our formalism, we use the extended Press- 



Figure 10. The bias of HII as a function of redshift in the model 
with fi = —3 source spectrum, but assuming an unbiased source 
population. 

Schechter theory to model the source clustering and adopt 
three different power-law type source spectra. We find that 
for UV dominated source spectra, the biases of the HII re- 
gions and the UV photons remain high during most of the 
reionization process. For hard source spectra, the HII regions 
tend to be ionized in a more homogeneous manner. The HII 
bias decays faster with time due to a longer photon mean 
free path and due to secondary ionization. The topology of 
the HII distribution is always inside-out, with overdense re- 
gions more highly ionized, at least on large scales. 

Our findings suggest that clustering measurements from 
future 21cm and CMB observations can be used to put con- 
straints on properties (both clustering and spectra) of the 
ionizing sources. Furthermore, on sufficiently large scales, 
both HII and HI have a scale independent bias with respect 
to dark matter - this means the same observations might be 
used to measure the shape of the matter power spectrum. 

A direct comparison of our results with 3D simulations 
(such as those by Kohler et al. 2005, Iliev et al. 2005, and 
Zahn et al. 2006) would be valuable and interesting. An 
accurate comparison requires detailed information on the 
source properties and the clumping factors from the simula- 
tions, and we hope to perform such a comparison in a future 
paper. 
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APPENDIX A - ANALYTIC SOLUTIONS TO 
THE FIRST ORDER RADIATIVE TRANSFER 
AND IONIZATION EQUATIONS 

Use eq.[23] in eq.[2D], we find the following equation for 



8lu 
where 



S eff -FA HII + / du'K(u>,cj')Ami{w) (38) 
sin[P(cj, u)')k] 



P(uu,uu')k 

X / djJ,B(LU, /j.) (k,(h)} R(lo' , n + lu — lu') 



(39) 



x exp[— / duu" B(uj" , fi + lu — lu' 1 



and 

du'K(u,w')8((j') (40) 

- oo 

f u , ,shi[P(u/,Lu')k] AT . ,. 
+ 4vr / dJ — „, v . 1 N(lu ) 

t> oc 

X / dflB(LO, (i)(/t(fl))A s (lil', (A + LU — Lu') 



x exp[— / duj" B{lo" , n + lu — lu")] 



To solve eq.[3S], we use an integrating factor 8(lu) which 
is defined as: 



d0_ 

d^jj 



(lu)F(lu) 



(41) 



and multiply both sides of eq. [38] by 0(lu): 



d(eA Hn ) =e $ eff + Q I dcu ' K (u;,Lu')A HII (Lu') (42) 



— $eff + / dLuK{LU,Lu) 



'duU 

or 

d{9A H ii 

6dLU 



Defining the function f(cu) 
d(6A HII ) 



(43) 



X fT 1 ^') / dLU 



„d{6A HII ) 



Olu>- 



6dLU 



(44) 



and changing the order of the integration, we can re-write 
eq.[43] as: 



f(uj) = 8 efs {uu) + I dLU T(iU,LU )f(LU ) 
where 



(45) 



T{lu,lu') — I du>" K (lu , lu" ) 



6(lu") 



(46) 



Eq.[45] is a standard Volterra integral equation of the sec- 
ond kind. Its solution can be easily generated by inverting 
a triangular matrix. Using f(uj), it is not hard to calculate 
Ah ii- One can then use eq.[23] to get A 7 or eq.[24] to get 
the monopole of A 7 . 



APPENDIX B - ON THE VALIDITY OF 
LINEAR PERTURBATION THEORY 

In this Appendix, we address the following question: is lin- 
ear perturbation theory justified on large scales even when 
highly nonlinear structures exist on small scales? This is 
a deep question that arises both in the present context of 
reionization and in the more familiar context of large scale 
structure formation theory. We will make no attempt to pro- 
vide a rigorous justification here. Intead, we will make some 
plausibility arguments, which are borrowed from the field 
of large scale structure (Peebles 1980). Ultimately, numer- 
ical simulations are needed to rigorously justify the use of 
perturbation theory on large scales. 

The fundamental equations are eq. [3] for ionization bal- 
ance and eq.[4] for radiative transfer. If the right hand sides 
of these equations were zero, these equations simply express 
conservation for nan and n 7 . If this were to hold true, we 
expect perturbation theory to work on large scales just like 
it is known to do for large scale structure - similar conserva- 
tion equations appear in large scale structure. Let us instead 
focus on the right hand sides which contain the novel aspects 
of the reionization problem. Morally one can think of eq.[3] 
as : 



drti 



dr 



ionization — recombination 



(47) 



= source 



sink 



(48) 



and one can think of eq. [4] as: 

dn-y 
~~dV 

The 'ionization', 'recombination' and 'sink' terms all contain 
quadratic combinations of the dynamical variables and are 
potentially what could cause the break down of perturbation 
theory. Let us divide space into regions where perturbations 
are large and where they are small. We will refer to these as 
the 'linear' and 'nonlinear' regions. Taking Fourier transform 
of the equation for ionization balance, we have 



+ 



3 dURII ik-x 

d x — e 

dr 



d x\ ionization — recombination le ! 



d x \ ionization — recombination ]e" 



(49) 



The linear regions are simple to deal with: we can linearize 
the 'ionization' and 'recombination' terms. The nonlinear 
regions potentially could give large contributions, but here 
we make use of a key insight: such regions are often in ion- 
ization equilibrum (i.e. 'ionization' roughly balances 'recom- 
bination') making the nonlinear contributions actually quite 
small. A very similar argument is used in the context of large 
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scale structure, where one invokes virial equilbirum as op- 
posed to ionization equilibrium to argue for the cancelation 
of potentially large terms. What are these nonlinear regions 
in our case? They could be HII bubbles at the beginning 
of reionization, or self-shielded Lyman-limit systems at the 
end of reionization, for instance. 

Setting the last term of eq.[49] to be zero, we have 

3 driHii ik-3 /,.„% 

a x — - — e = (50) 
dr 

d 3 xW(x) [ ionization — recombination ]ii ncar e 8fe :E 

where we have introduced a mask W(x) which vanishes in 
the nonlinear regions, and equals unity in the linear regions, 
and we have linearized the ionization and recombination 
terms. 

Ultimately, we are interested in the power spectrum of 
fluctuations in uhii for instance. Eq. (|50[) tells us that its 
evolution can be regarded as linear as long as we are study- 
ing scales k on which the mask W, or its Fourier trans- 
form, has a negligible effect on the power spectrum. Gener- 
ally, the linear approximations will be valid only on scales 
much larger than the size of the nonlinear regions, and likely 
even larger than the clustering scale of the nonlinear regions. 
For instance, in the early stages of reionization, we expect 
perturbation theory to work only on scales that encompass 
many HII bubbles. Note that the mask in question is more 
complex than the usual masks in galaxy surveys: here, the 
mask is correlated with the signal in non-trivial ways, and it 
is by no means obvious that the large scale power spectrum 
is unaffected by such masking. Addressing this important 
issue is beyond the scope of this paper. 

How about the radiative transfer equation? Applying a 
similar split, we have 



3 dn^ ik-x /ri \ 

d x— — -e = (51) 
dr 



d 3 :r[source — sinkje* 



+ / d ^[source — sinkje 1 " x 

J nonlinear 

Here, the situation is similar to the ionization balance equa- 
tion, in that in nonlinear regions, one expect the 'source' 
and 'sink' to roughly cancel, but in general, we don't ex- 
pect exact cancelation. For instance, one can think of a 
galaxy where the UV photons are propagating out of a 
thick medium. Most of the photons are consumed within 
the galaxy, but inevitably there will be some that escape. 
One usually quantifies this by the escape fraction. Here, we 
can account for this effect by renormalizing the source: 

fodnyji.a = (52) 

dr 

d 3 a;[ source rcnorm. — W(x) sinkje* ' x 

The mask W is the same as before, and the 'sink' term can 
be linearized. The 'source' term can be thought of as a new 
effective source that accounts for the escape fraction from 
the nonlinear regions. 
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